/*
 
 Trabalho 1
 Programação Concorrente
 
 Cálculo de PI com precisão de 10 milhões de dígitos
 Método: Borwein (Paralelo)
 Biblioteca de precisão: The GNU Multiple Precision Arithmetic Library (GMP)
 Biblioteca para threads: POSIX thread functions
 
 */

#include <stdio.h>
#include <gmp.h>
#include <pthread.h>

#define LOG_2_10 3.32192809488736234787

/* Variáveis big float da GMP */
mpf_t pi_old, pi_new;
mpf_t a_old, a_new;
mpf_t y_old, y_new;
mpf_t b, c, d, e, f, g;
mpf_t epsilon;
/* Fator multiplicativo */
int k = 8;

/*
 
 As funções abaixo rodam em uma thread separada cada
 para calcular partes independentes da fórmula
 
 */

void * calculate_c (void * args)
{
    mpf_ui_sub(c, 1, b);
    return NULL;
}

void * calculate_d (void * args)
{
    mpf_add_ui(d, b, 1);
    return NULL;
}

void * calculate_f (void * args)
{
    mpf_pow_ui(f, e, 4);
    mpf_mul(f, f, a_old);
    return NULL;
}

void * calculate_g (void * args)
{
    mpf_pow_ui(g, y_new, 2);
    mpf_add(g, g, e);
    mpf_mul(g, g, y_new);
    mpf_mul_ui(g, g, k);
    return NULL;
}

int main (void)
{
    /* Contador de iterações */
    int i = 1;
    /* variável para gravar o nome do arquivo de saída */
    char filename[10];
    /* ponteiro para arquivo de saída */
    FILE * F;
    /* declaração das threads */
    pthread_t thread_c, thread_d, thread_f, thread_g;
    
    /*
     
     Define a precisão padrão usada nos cálculos:
     
     OBS: Utilizamos 10000100 casas de precisão para
     possíveis overflows durante o cálculo.
     Como o argumento é em bits, e não em casas decimais,
     multiplicamos por uma constante log2(10),
     que é a quantidade média de bits por casa decimal
     
     */
    mpf_set_default_prec(LOG_2_10 * 10000100);
    
    /*
     
     Inicializa as variáveis de acordo com o método:
     a = 6-sqrt(32)
     y = sqrt(2)-1
     pi = 1/a
     
     */
    mpf_init(pi_old);
    mpf_init(pi_new);
    mpf_init(a_old);
    mpf_init(y_old);
    mpf_init(b);
    mpf_init(c);
    mpf_init(d);
    mpf_init(e);
    mpf_init(f);
    mpf_init(g);
    
    mpf_init_set_d(a_new, 32.0);
    mpf_sqrt(a_new, a_new);
    mpf_ui_sub(a_new, 6, a_new);
    
    mpf_init_set_d(y_new, 2.0);
    mpf_sqrt(y_new, y_new);
    mpf_sub_ui(y_new, y_new, 1);
    
    mpf_init_set_str(epsilon, "1e-10000000", 10);
    
    mpf_ui_div(pi_new, 1, a_new);
    
    /*
     
     Calcula as demais iterações, gravando cada resultado parcial
     em um arquivo de nome it-x.txt
     
     */
    do
    {
        mpf_swap(pi_old, pi_new);
        mpf_swap(a_old, a_new);
        mpf_swap(y_old, y_new);
        
        mpf_pow_ui(b, y_old, 4);
        mpf_ui_sub(b, 1, b);
        mpf_sqrt(b, b);
        mpf_sqrt(b, b);
        pthread_create(&thread_c, NULL, calculate_c, NULL);
        pthread_create(&thread_d, NULL, calculate_d, NULL);
        pthread_join(thread_c, NULL);
        pthread_join(thread_d, NULL);
        
        mpf_div(y_new, c, d);
        mpf_add_ui(e, y_new, 1);
        
        pthread_create(&thread_f, NULL, calculate_f, NULL);
        pthread_create(&thread_g, NULL, calculate_g, NULL);
        pthread_join(thread_f, NULL);
        pthread_join(thread_g, NULL);
        k *= 4;
        mpf_sub(a_new, f, g);
                
        mpf_ui_div(pi_new, 1, a_new);
        
        mpf_sub(pi_old, pi_new, pi_old);
        mpf_abs(pi_old, pi_old);
        sprintf(filename, "it-%d.txt", i++);
        F = fopen(filename, "w");
        gmp_fprintf(F, "%.10000000Ff\n", pi_new);
        fclose(F);
    } while ( mpf_cmp(epsilon, pi_old) < 0 );
    
    /* Libera o espaço de memória alocado para as estruturas */
    mpf_clear(pi_new);
    mpf_clear(pi_old);
    mpf_clear(a_new);
    mpf_clear(a_old);
    mpf_clear(y_new);
    mpf_clear(y_old);
    mpf_clear(b);
    mpf_clear(c);
    mpf_clear(d);
    mpf_clear(e);
    mpf_clear(f);
    mpf_clear(g);
    mpf_clear(epsilon);
    
    return 0;
}
